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Abstract 

We present a numerical study of a reaction-diffusion model on a small-world network. We 
characterize the model’s average activity Ft after T time steps and the transition from 
a collective (global) extinct state to an active state in parameter space. We provide an 
explicit relation between the parameters of our model at the frontier between these states. 
A collective active state can be associated to a global epidemic spread, or to a persistent 
neuronal activity. We found that Ft does not depends on disorder in the network if 
the transmission rate r or the average coordination number K are large enough. The 
collective extinct-active transition can be induced by changing two parameters associated 
to the network: K and the disorder parameter p (which controls the variance of A). We 
can also induce the transition by changing r, which controls the threshold size in the 
dynamics. In order to operate at the transition the parameters of the model must satisfy 
the relation rK = Op, where ap as a function of p/{I — p) is a stretched exponential 
function. Our results are relevant for systems that operate at the transition in order 
to increase its dynamic range and/or to operate under optimal information-processing 
conditions. We discuss how glassy behaviour appears within our model. 
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Many problems in Science can be cast in terms of dynamics on networks: social 


phenomena [l|j_|2 
brain activity 


idemic spread [3|, food we bs B and ecosystem’s diversity 


11 , 


H, granuW materials UMM and, in general, complex 


systems 0. Among the most studied models, the small-world network model of Watts 
and Strogatz (WS) [13, EM 0 can be tuned to interpolate between a regular and a 
random network, a very atractive property that allows us to explore the consequences of 
network disorder on dynamics. In this work we consider a stochastic reaction-difusion 
cellular automata model on a small-world network and study its average activity after 
T time steps and its collective extinct-active transition. For the first time, we provide a 
explicit relation between the parameters of the model for the system to operate at the 
transition, and disorder enters in this relation as a stretched exponential function. 
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In the WS model we start from an ordered ring in which each node is conected to its 
K nearest neighbours, and then we randomly rewire each conexion with probability p. 
Since we only rewire and we don’t create new conexions, the average coordination number 
is still K, but the distribution of K gets broader as we increase p 18|. Two quantities 
typically used to characterize a network are the mean minimal distance between two 
nodes L and the clustering coeficient C. For each node, we can measure how many of it’s 
neighbours are conected between them. C is defined as the actual number of conexions 
between a node’s neighbours over the maximun number of posible conexions between 
neighbours, averaged over all nodes. Hence, in a social network C could measure up to 
what extent someone friends are friends between them. In the WS network there is a 
range of values of p for which C is still high and L is small, a signature of a small-world 
scenario. The case p = 0 corresponds to an ordered network, with hi gh C and L. The 


case p = 1 corresponds to a random network, with small C and L 18|, Il7 |. 


A very generic reaction-difusion model is the Greenberg-Hastings model (GH). This 
model has been used to model the Belousov chemical reaction, biolo gica l systems, epi- 
demic spread, collective behaviour of heart cells and neuronal activity 1^, 20,[2ll[22, 23, 
24 12^ It]. Our stochastic version of this three state model is as follows: 


• If a cell is in the excited state at time t then it is in the passive state at time t + 1. 

• If a cell is in the passive state at time t then it is in the susceptible state at time 
t 1. 

• If a given cell is in the susceptible state at time t, and at least one of its neighbours 
is in the excited state at time t, then the given cell is in the excited state at time 
t -f 1 with probability r, otherwise the given cell remains in the susceptible state. 

The GH is a model for excitable media ([1^, see fig. 1 in [^) whose threshold is 
controlled by the parameter r. We call the parameter r the transmission (or infection) 
probability. If excited means no-healthy, then I become infected if at least one of my 
neighbours is infected, but with probability r. In this context r would be a measure 
of the average state of the inmune system of the population, with a smaller r implying 
more resistance to become ill. If we are in a social context, r could be associated to 
the confidence of the social agents on an especific action, with greater r implying more 
confidence. For neuronal activity or chemical reactions, r can be associated to a potential 
barrier, with a smaller r implying a greater potencial barrier (a greater threshold). In 
a raw model for dense granular flows, the states of the GH model could be associated 
to fluid or solid type of contacts between grains 27|, with the excited state associated 
to fluid contacts, the passive state associated to contacts with static friction and with a 
large time of contact 271 or with a small mobilization of friction [l^, and the susceptible 
state associated to contacts with static friction and with a short time of contact or with 
a large mobilization of friction. Versions of the GH model have been implemented on 
regular networks 28|, [22|, [24 1, scale free networks 1^ , on a small-world network to study 
spiral waves and on the human connectome Q- 

It is of particular interest under what conditions the system, for a given initial condi¬ 
tion, evolves towards a global active state, in which a finite fraction of the nodes remains 
active 22|. A collective active state can be associated to a global epidemic spread, or to 


a persistent neuronal activity. For random initial conditions, we have found transitions 
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Figure 1: Average activity Ft as a function of the average coordination number A, for several values of 
p and r. In the inset is shown the plateau of Ft-, obtained for large A, as a function of 1/r. For r = 0.05 
and p = 0.1, for example, we have a collective active state for K > 32. If r ^ 1 or if A is large enough 
Ft becomes independent of the disorder parameter p. The number of nodes is N = 1000 and we show 
the average result for 100 realizations for each set of parameters. For this figure, we looked for activity 
in the system after T = 1000 time steps. 
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to a collective active state as we vary any of the three parameteres in the GHWS model: 
the transition can be induced by increasing the average coordination number K, by in¬ 
creasing the disorder in the network p or by increasing the transmission probability r. 
Based on our numerical results, we provide for the first time a explicit relation between 
the parameters of the model for the system to operate at the transition, and disorder en¬ 
ters this relation as a stretched exponential function. Stretched exponentials are usually 
found in glas ses and disordered systems [iillllsi. As far as we know, for this model 
(GHWS) [291 we are presenting a novel way to consider the influence of disorder on the 
extinct-active frontier. In reference 0 the network was the human connectome and only 
a threshold was varied in order to tune the transition. 

Our system consists of = 1000 nodes. Initially, each node is randomly asigned, 
with equal probability, to one of the three states of the model. We call F the average 
number of nodes in the excited state. The behaviour of the GHWS model is sumarized 
in figure [TJ on which we show the average activity after T = 1000 time steps Ft{K,p, r). 
We have a collective active state if after T time steps (see below) Ft > 0. We have a 
collective extinct state if after T time steps Ft = 0. For K > Kc we have a collective 
active state, and from this figure we see that Kc = Kc{r,p). It can be observed that 
for r —>■ 1 Ft becomes independent of the disorder parameter p. Also, for any given 
transmission probability r, Ft becomes independent of p if AT is large enough: it tends 
to a plateau whose dependence on 1/r is shown in the inset. 

For r = 0.05, we show in figure [5^ zones of global active states in parameter space. 
We have found, as previously reported [s^, that for K < Kc the relaxation time of decay 
to zero activity increases largely as we approach a transition to a collective active state. 
Because of this, we looked for activity after T time steps, and by a time step we mean an 
update of the whole network. Thus, we are characterizing the evolution in time of the 
frontier, and the limit T ^ oo corresponds to the usual meaning of the frontier. For a 
given r, we can induce a transition to a collective active state by increasing the disorder 
parameter p or by increasing the average coordination number K. 

In order to characterize the frontier between active and extinct global states in pa¬ 
rameter space we show in figures [5 ]d and [5J: the critical transmission probability as a 
function oil/K, for disorder spanning four orders of magnitude in p. Our results can be 
sumarized in the following simple relation: 


K' 


( 1 ) 


This scaling of Tc with K was obtained by Berry and Fates [2^ in a mean field 
approximation of a similar GH model, where the crucial ingredient was to obtain an 
approximation to the conditional probability that in the neighbourhood of a given node 
there is at least one node in the excited state, given that the considered node is in the 
susceptible state. The scaling Vc ^ K~^ was verified by numerical simulations on a 
regular network by considering different types of neighbourhoods, obtaining Cp « 2 (the 
exact mean field result was Op = I) [1^; this result was robust against the inclusion of 
defects in the regular network. In figure[2ji we show Op as a function of p, and we conclude 
from this figure that disorder in the network reduces Tc- Since we have a collective active 
state for r > within the GHWS model disorder favours a collective active state. 

In figure [3] we see a detailed version of figure [3Ji. Our numerical results for Op are 
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Figure 2: a) Zones of collective active states for r = 0.05 after T time steps, b) Critical transmission 
probability Tc as a function of 1/A", estimated after T = 100 time steps for different p. c) Same as in 
b) but with T = 1000. As in b), the dashed lines corresponds to the relations Vc = l/A" and Vc = 2/K 
. Our system is a network with N = 1000 nodes and we averaged over m = 100 realizations. Dotted 
lines, here and in b), are best fits to the relation Vc = apiK. We don’t considered the case r = 1 (the 
deterministic GH model) for which we have a collective active state for A" = 2, for any p. d) The slope 
ap as a function of the disorder parameter p, obtained from c) and b). 
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Figure 3: This figure is a detailed version of figure [2]i. In order to obtain each point in this graph we 
considered for each p five values of K\ iC = 4, 8,16,32,64, with m realizations for each K. We then 
adjusted the best line to equation Q. For T = 50 and T = 100, we used m = 1000; for T = 1000 we 
used m = 100; for T = 10000 we used m = 30. Dashed lines are best fits to stretched exponentials of 
the form given by equation (2] In the inset we show the four parameters of eq. [2] as a function of T. The 
dot-dashed line is equation [J] for T = 1000. We show, in stars, the results for N = 10000 and T = 100, 
with no major differences with the case N = 1000; certainly L depends on N [T3l but L appears to affect 
only the time required to propagate the activity to the whole system El- The dotted line is the mean 
field result ap = 1 [2^ . Recently (v6), I have been aware that if we average in a different way we do have 
a N dependence: if we do only one realization for each value of K (see the beginning of this caption), 
obtaining ap and then average several realizations we do have a N dependende in the results. This fact 
will be considered elsewhere. For N = 1000, the results presented here do not change qualitatively. 
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Figure 4: a) Behaviour of a; = (up — ai)/A, and it’s derivative, with x = ^dw/dz. b) ap for T = 10®, 
T = 10^ and T = 10® (different run from previous figure). Dashed lines are best fits to stretched 
exponentials of the form given by eq. l(2]|. At T 10^ the frontier no longer changes it’s shape, all the 
parameters of the stretched exponential no longer changes with time, except for ai. For rK < ap we 
have an extinct state, c) The frontier keeps moving since ai keeps increasing (as ai 0.015 ln(T) for 
T > 2000). We could say that for T > 2000 the frontier enters a glassy-like regime. 


consistent with a stretched exponential of the form: 

dp = Aexp[-(a:/T)^] + ai (2) 


with X = p/{1—p). Stretched exponentials are usually found in relaxation of glasses and 
disordered systems [s^ [sl, [s^ [sB, ll^l • In the inset of figure [3] we show the evolution of 
the four parameters of equation [2] 

In terms of z = ln(l/a:), there is a maximum in the quantity x(-®) = dapjdz at 
Zc = ln(l/r), with the maximum given by Xm = xi^c)' 



(3) 


The quantity x measures the sensitivity of the extinct-active frontier to changes in 
disorder (through the variable z), going to zero for very ordered or very disordered 
networks. In terms of /i = e~^ jx, the order parameter w = (op — ai)/A = exp(— 
and x(z) = IS^dwIdz, we have: x{^)IXm = e/iic (see figure 0^). If we operate at Zc we 
can communicate more easily through collective states when considering, for example, 
a network of networks 38| in which each unit is a GHWS one: at the frontier and at 
Zc small changes in disorder can allow the system to switch between an extinct and an 
active collective state more easily. 

If we expand around Zc, we obtain the approximate relation: 


Op - Qpc ^ Xm(z - Zc) (4) 

with Opc = A/e + ai- In figure [3] we show eq. 0]for T = 1000. 

By monitoring the activity T’ as a function of T we have found that Op is well described 
by a stretched exponential function in the variable x = p/(l — p), for all T. Since at 
the extinct-active frontier rK = Up (eq. [T]), we can interpret Op as the average minimal 
number of effective neighbours to become active. In eq. [2] we have defined A as the 
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difference between ag (up when p —>■ 0) and ai (up when p 1), see the inset of figure [S] 
What happens, for example, if we consider an ensamble of systems such that rIC ss 1.63, 
in which each member of the ensamble can have any value of p?. We see in figure |3] that 
oo ~ 1.63 for T = 100. Thus, up to T = 100 we’ll have ordered and disordered members 
of the ensamble that are still active, but in the long run only the more disordered members 
of the ensamble (x > 0.2, approx., see fig. ED will stay active. If we consider an ensemble 
of systems with rK < oi, then in the long run all the members of the ensemble will go 
extinct. On the other hand, if we consider an ensemble of systems with rK > ao then 
all the members of the ensemble will stay active in the long run. 

In figure |3 )d we can see that the frontier no longer changes it’s shape for T > 10^, 
only oi keeps increasing with time, as can be seen in figure 0J:. With the relation oi = 
1.18 + 0.015 ln(T) being satisfied for T > 2000 (see fig. IH:), we have that for an ensamble 
with rK Ki 2.4, at Tg = 10® the more ordered members of the ensamble starts to go extinct 
(see fig. Sb). For this ensamble, if we denote Ti as the time when the more disordered 
members of the ensamble will go extinct we have that TijTg = k, ~ 10^®. 

This is reminiscent of glassy behaviour. 

Several generalizations can be introduced in the model implemented here. One of 
them would be to consider a weighted network, with weights that can be correlated or 
not to the local coordination number Q. In order to consider different passive time 
scales, a possible generalization of our GH model would be to introduce M time steps in 
the passive (refractory) state before becoming susceptible. Regarding this possibility, in 
reference [22| was found that, for a regular network, the critical transmission probability 
depends weakly on M, and that at the active-extinct transition the decay in time of the 
average activity F follows a power law, with an exponent that does not depends on M: 
its universality class remains directed percolation. 

We have characterized the activity of a GHWS model, and we found that the average 
activity does not depends on disorder in the network if the transmission rate r or the 
average coordination number K are large enough. We have found that a collective extinct- 
active transition in a GHWS model can be induced by increasing r, a parameter that 
controls the threshold size in the system’s dynamics. Also, the transition can be induced 
by increasing the disorder parameter p and by increasing K. Our results are relevant 
to systems that need to operate at the extinct-active transition, in order to increse 
its dynamic range and/or to operate under optimal information-processing conditions 
(3ll| . Remarkably, in reference Q it was found that in order to reproduce the patterns 
of neuronal activity observed in experiments, the threshold of a stochastic GH model 
implemented on the human connectome has to be tuned for the system to operate at 
the the extinct-active transition [s^. Within the GHWS model studied here, in order to 
operate at the transition the parameters of the model must satisfy the relation rK = ap, 
where ap as a function oi p/{I — p) is a stretched exponential function. 

I would like to acknowledge Miguel Pineda for usefull discussions during the first 
stages of this investigation. 
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